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2-D Radiative Transfer in Protostellar Envelopes: I. Effects of 

Geometry on Class I Sources 
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ABSTRACT 

We present 2-D radiation transfer models of Class I Protostars and show the 
effect of including more realistic geometries on the resulting spectral energy dis- 
tributions and images. We begin with a rotationally flattened infalling envelope 
as our comparison model, and add a flared disk and bipolar cavity. The disk 
affects the spectral energy distribution most strongly at edge-on inclinations, 
causing a broad dip at about 10 /iin (independent of the silicate feature) due to 
high extinction and low scattering albedo in this wavelength region. The bipolar 
cavities allow more direct stellar+disk radiation to emerge into polar directions, 
and more scattering radiation to emerge into all directions. The wavelength- 
integrated flux, often interpreted as luminosity, varies with viewing angle, with 
pole-on viewing angles seeing 2-4 times as much flux as edge-on, depending on 
geometry. Thus, observational estimates of luminosity should take into account 
the inclination of a source. The envelopes with cavities are significantly bluer 
in near-IR and mid-IR color-color plots than those without cavities. Using 1-D 
models to interpret Class I sources with bipolar cavities would lead to an un- 
derestimate of envelope mass and an overestimate of the implied evolutionary 
state. We compute images at near-, mid-, and far-IR wavelengths. We find that 
the mid-IR colors and images are sensitive to scattering albedo, and that the 
flared disk shadows the midplane on large size scales at all wavelengths plotted. 
Finally, our models produce polarization spectra which can be used to diagnose 
dust properties, such as albedo variations due to grain growth. Our results of 
polarization across the 3.1 /im ice feature agree well with observations for ice 
mantles covering 5% of the radius of the grains.. 
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1. Introduction 

Lada and collaborators first classified pre-main sequence stars into an evolutionary se- 
quence, Class 1 through III, by the shapes of their spectra energy distributions (SEDs) (Lada 
1987; Myers et al. 1987; Adams, Lada, & Shu 1987). Class I sources, the subject of this 
paper, are thought to be T Tauri stars surrounded by infalling envelopes. More recently, 
an earlier stage was proposed, the Class stage to describe more deeply embedded sources 
(Andre, Ward-Thompson, & Barsony 1993). This called into question whether the Class I 
envelopes are merely remnant or are still collapsing, though detailed radiative transfer mod- 
els of several Class I sources require substantial envelope mass (typically 10-20 times greater 
than a Class II disk), and suggest that the envelopes are still infalling, if at a less embedded 
stage than the Class stage (Adams et al. 1987; Kenyon et al. 1993a,b; Whitney, Kenyon, 
& Gomez 1997; Stark et al. 2002). 

Adams et al. (1987) successfully interpreted the SEDs of Class I sources with 1-D 
radiative transfer models of spherically-averaged rotationally-flattened infalling envelopes. 
Kenyon et al. (1993a) applied "1.5-D" radiative transfer calculations to model the SEDs of 
the Taurus Class I sources, using the spherically-averaged density for calculating the tem- 
perature. The average model infall rates and centrifugal radii were consistent with observed 
sound speeds and disk sizes. Efstathiou & Rowan-Robinson (1991) also performed accurate 
2-D radiative transfer solutions using the rotationally-flattened infalling envelope. All of 
these SED models of Class I sources tend to have similar features, steeply rising at about 1 
/xm, then gradually rising to 100 /xm, and falling off to 1000 /im (see Figure 9, top left panel 
for a similar model). More recently, Nakazato, Nakamoto, & Umemura (2003) calculated 
2-D models of Class I objects which included a disk and a piece-wise spherically symmetric 
envelope divided into two latitude sectors, each having a radial power law independent of 
polar angle. To create a bipolar cavity, the density was lowered in the polar sector. 

Most of the previous 1-D and even the 2-D SED models of Class I sources underestimate 
the optical/near-IR flux when compared to observations. Scattered light models (Kenyon 
et al. 1993b) showed that bipolar cavities are required to allow an escape path for the 
short-wavelength radiation, and to fit the extended images and polarization maps of these 
sources (Whitney et al. 1997, Lucas & Roche 1997, 1998). Gomez, Whitney, & Kenyon 
(1997) found [S II] emission from jets and Herbig-Haro objects in nearly all the Taurus Class 
I sources, providing a mechanism for producing the bipolar cavities. In addition, most of 
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the Class I sources are known to have low-velocity molecular outflows (Moriarty-Schieven 
et al. 1992, 1995a,b; Chandler et al. 1996, Tamura et al. 1996, Hogerheijde et al. 1998). 
High spatial-resolution images (Padgett et al. 1999) and scattered-light models (Stark et 
al. 2002) find that vertically extended disks, an expected consequence of rotational infall, 
are required in addition to the infalling envelopes to fit the images several Class I sources. 
The hydrodynamic collapse models of Bate (1998) also produce substantial disks in Class I 
sources. 

Since it is clear from these studies that Class I sources have disks and cavities in addi- 
tion to infalling envelopes, we wish to determine if these more complicated geometries give 
different SEDs than the previous models show. This paper computes the exact solutions for 
2-D radiative transfer in Class I sources including the realistic geometrical effects of a rota- 
tionally flattened envelope, a flared accretion disk, and bipolar cavity. We find that both the 
disks and cavities have substantial effects on the resulting SEDs, especially in the mid-IR. 
These variations can be tested with new ground-based and upcoming space-based infrared 
observations. In addition to SEDs, we compute images and polarization. This paper details 
how the geometry affects the SED for a single evolutionary state, the Class I source. In a 
followup paper we will explore these effects in an evolutionary sequence. The next section 
describes the models, §3 shows the results of four different geometries, and §4 concludes with 
a discussion of the results. 



2. Radiative Transfer Models 

2.1. Ingredients of the Model: Disk, Envelope, & Cavity Structure 

We wish to use as realistic a geometry as possible to describe the Class I protostar. 
Since disks appear to be ubiquitous in star forming regions (Cohen, Emerson, & Beichman 
1989; Strom et al. 1989; Beckwith et al. 1990, Dutrey, Guilloteau, & Simon 1994; Koerner 
& Sargent 1995; Saito et al. 1995; Burrows et al. 1996), and therefore the infalling envelope 
likely has angular momentum, we choose for the envelope geometry the infall solution that 
includes rotation (Ulrich, 1976, Cassen & Moosman 1981, Terebey, Shu, & Cassen 1984). On 
larger scales, magnetic fields (Galli & Shu 1993a,b, Li 1998, Ciolek & Konigl 1998) and initial 
motions prior to dynamical collapse (Foster & Chevalier 1993) may affect the structure but 
within the infalling region, the effects of rotation likely dominate over magnetic effects. This 
region contributes most to the SED and images. The density structure for the rotationally- 
flattened infalling envelope within the infalling radius is given by: 
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where M env is the envelope mass infall rate, R c is the centrifugal radius, /i = cos 6 and fi is 
the cosine polar angle of a streamline of infalling particles as r — > oo). The equation for the 
streamline is given by 

l4 + Mr/Rc - 1) - ^(r/i?c) = . (2) 

For the disk structure, we use a standard flared accretion density (Shakura & Sunyaev 1973, 
Lynden-Bell & Pringle 1974; Pringle 1981; Bjorkman 1997; Hartmann 1998): 
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where w is the radial coordinate in the disk midplane and the scale height increases with 
radius, h = h [vo / R*) 13 . For the models presented here, we adopt: (3 = 1.25, a = 2.25, 
ho = O.Oli?*, giving a scale height at the disk outer radius of /i[100AU] = 10AU. The f3 
and a values are consistent with detailed hydrostatic structure calculations (D'Alessio et 
al. 1998). The exact disk structure in an object with an infalling envelope is likely to 
vary from this geometry, since the disks are accreting from the infalling envelope. However, 
observations and models of Class I sources show morphologies consistent with vertically 
extended disks (Padgett et al. 1999, Stark et al. 2002) so the flared disk seems to be a 
good starting point for this investigation. We include accretion energy, dE^/dV, in the 
disk using the a— disk prescription (Shakura & Sunyaev 1973, Kenyon & Hartmann 1987, 
Bjorkman 1997; Hartmann 1998): 
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where the disk accretion rate, M disk is given by 

M disk = a disk V c po hl/R* , (5) 

with the critical velocity V c = ^/GM^/R*. For the disk parameters used in this paper (§3), 
Table 1 shows the resulting disk accretion rate and accretion luminosity. For simplicity, ac- 
cretion energy inside the dust destruction radius, Rdust, is emitted with the stellar spectrum. 
Thus the accretion luminosity emitted within the dusty disk is 
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Since near-IR images of Class I sources require bipolar cavities to fit the observations 
(Kenyon et al. 1993b, Whitney et al. 1997), we include them in these models. We choose two 
shapes to describe the cavities since they both appear to be evident in observations (Padgett 
et al. 1999, Stark et al. 2002; Reipurth et al. 2000). These are a streamline, which is conical 
on large scales, and a curved cavity. The streamline cavity shape might occur if precessing 
jets carve out a conical shape while infalling material outside of the cavity continues to fall in 
along streamlines. The cause of a curved shape is less certain. The streamline cavity shape 
is defined by equation (2); at each (/x, r), /i is calculated, and if it is greater than the chosen 
cavity angle, the density is set to the cavity density. Whitney & Hartmann (1993) show a 
picture of this cavity shape in their Figure 1. The curved cavity follows z = aw b , where 
zu = x 2 + y 2 . We include a small amount of dust in the cavity with constant density, 
nu2 = 2 x 10 4 cm~ 3 , as expected for a cylindrical outflow. The density in the cavity may 
vary from nui ~ 10 3 — 10 6 cm -3 in the high- and low- velocity flow regions (Bachiller & 
Tafalla 1999; Reipurth & Raga 1999). A more realistic cavity density distribution including 
such large variations could have a substantial effect on the radiative transfer, and will be 
considered in future work. The total mass in the cavity, about 6 x 1O" 3 M , is similar to 
observed outflow masses in Taurus Class I sources (Moriary-Schieven et al. 1995b, Chandler 
et al. 1996, Hogerheijde et al. 1998). 

2.2. Radiation Transfer Method 

To calculate the radiative transfer, we incorporated the Monte Carlo radiative equilib- 
rium routines developed by Bjorkman & Wood (2001) into a 3-D spherical-polar grid code 
(Whitney & Wolff 2001). Bjorkman & Wood tested their codes by comparing to a set of 
benchmark calculations developed by Ivezic et al. (1997) for spherically symmetric codes. 
They then applied their method to 2-D ellipsoidal envelopes and T Tauri disks (Wood et 
al. 2002a,b). These codes use a numerical optical depth integrator that operates on an an- 
alytically prescribed density. The numerical integrator has difficulty handling sharp density 
gradients. In order to include arbitrary density distributions, such as bipolar cavities or 3-D 
clumpy envelopes, we have developed a code that uses a 3-D spherical-polar grid (Whit- 
ney & Wolff 2001). The density is constant within a grid cell, allowing for straightforward 
integration of optical depth through the grid. Since the density in each cell is computed 
only once at the beginning of the processing, we can incorporate very complicated analytic 
density formulations, or use tabular densities derived from hydrodynamical simulations. For 
the protostellar envelope case, the grid has variable spacing in r and 9 to sample the vastly 
different density variations and size scales from the inner disk region to the outer envelope. 
The radial spacing is logarithmic in the inner region of the disk, to resolve the inner edge 
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of the disk in optical depth, and then follows a power-law spacing out to the edge of the 
envelope, chosen to be 5000 AU for the models presented here. 

We include accretion in our disk, following Wood et al. (2002a). We also attach a 
diffusion solution to the inner regions (Bjorkman et al. 2003). We tested our grid code 
by comparing calculations of ellipsoidal envelopes and T Tauri disks to the 2-D codes of 
Bjorkman & Wood (2001) and Wood et al. (2002a,b). The results were identical. 

In these models, the luminosity is produced by the central star and disk accretion; it 
then gets absorbed and reemitted, or scattered, by the surrounding disk and envelope. As 
described in Bjorkman & Wood (2001), the radiative equilibrium algorithm corrects for the 
temperature in each grid cell by sampling the new photon frequency from a differenced 
emissivity function (the difference between that computed with the previous erroneous tem- 
perature and that computed from the new corrected temperature). Energy is conserved and 
no iteration is required as the difference spectrum corrects for the previous temperature er- 
rors. The scattering phase function is approximated using the Henyey-Greenstein function, 
with forward-scattering parameter, g, and albedo, oj computed at each wavelength from our 
grain model (described below). Polarization is calculated assuming a Rayleigh-like phase 
function for the linear polarization (White 1979). 

The Monte Carlo codes produce results simultaneously at all inclinations. We also incor- 
porate the "peeling-off" algorithm (Yusef-Zadeh, Morris, & White 1984, Wood & Reynolds 
1999) to produce high S/N images at a specific inclination in a given run. Images are pro- 
duced at several wavelengths and convolved with broadband filter functions for comparison 
to observations. For this paper, we chose for the near-IR images the Hubble Space Telescope 
NIC2 F110W (1.1 /im), F160W (1.6 /im), and F205W (2.1 /im) filters. For the mid-IR 
images, we chose the SIRTF IRAC filter functions at 3.6, 4.5, 5.8, and 8.0 /im. And for the 
far-IR images, we chose the SIRTF MIPS filter functions at 24, 70, 160 /im. A given run 
produces all the images simultaneously as it produces the SEDs and polarization spectra. 
The models take ~ 30 minutes to produce high quality spectra (that is, high frequency res- 
olution and low noise) at all inclinations on a 1 GHz linux PC running g77 FORTRAN. To 
produce high quality images takes several hours. 



2.3. Dust Model 

The dominant opacity source in the radiative transfer models is dust grains. Our dust 
model is based on that derived by Kim, Martin, and Hendry (1994) for the diffuse ISM 
(using a code kindly supplied by Peter Martin and Paul Hendry). To approximate dust 
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properties in nearby star formation regions, we produce a size distribution of grains that 
fits an extinction curve generated by the Cardelli, Clayton, & Mathis (1989) prescription, 
with Ry=4. This value of Ry is typical of the denser regions in the Taurus molecular 
cloud (Whittet et al. 2001). The size distribution is not a simple power law, because it is 
derived from a Maximum Entropy Method solution. The size distribution is similar to that 
shown in Figure la of Kim et al. (1994). We adopt the dielectric functions of astronomical 
silicate and graphite from Laor & Draine (1993). To illustrate the effects of ices, which 
cause observable spectral features, we include a layer of water ice on the grains, covering the 
outer 5% of radius. We varied this value between and 15% and found 5% to agree best 
with polarization observations, as discussed in §3.5. The ice dielectric function is taken from 
Warren (1984). 

This grain model may be appropriate for the envelope properties but is likely not appli- 
cable in the midplane of the disk where grains should grow to meter-size within a timescale 
of 10 4 - 10 5 years (Weidenschilling & Cuzzi 1993; Wetherill & Stewart 1993; Beckwith, Hen- 
ning, & Nakagawa 2000). Thus, at very long wavelengths where thermal emission from the 
disk dominates, our SEDs may not agree with observations (see D'Alessio, Calvet, & Hart- 
mann 2001, Wood et al. 2002a for models of disks using large grains). A plot of the grain 
properties vs. wavelength is shown in Figure 1. The properties for ISM grains are also shown 
for comparison. The major difference is that the albedo of the i?y=4 grains is substantially 
higher out to about 10 pm, which will give brighter images at these wavelengths than ISM 
grains. 



3. Four Variations on a Class I Source Geometry 

To illustrate the effects of geometry on the SED we will show results for four geometries. 
The first is the rotationally-flattened infalling envelope only, illuminated by the central star. 
This is the geometry used by several previous authors (though spherically averaged in some 
cases). For model parameters, we use typical values derived from previous models of Taurus 
sources (Kenyon et al. 1993a,b): infall rate, M = 5 x 10~ 6 M & /yr , centrifugal radius 
R c = 100 AU, envelope outer radius R ma x = 5000 AU, and envelope inner radius equal to 
the dust destruction radius, about 7.5 stellar radii in this case. The central star radiates as a 
4000 K star, using a Kurucz model atmosphere (Kurucz 1994), with a radius of R± = 2.O9i? , 
giving a stellar luminosity of 1 L Q . An azimuthal cut through the envelope density structure 
is shown in the top-left panel of Figure 2. It has a toroidal shape in the mid-plane about 
twice the size of the centrifugal radius. The density profile within this region falls roughly 
as p oc r -1 / 2 (see equation 1), and beyond this changes to p oc r -3 / 2 . 
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For our second geometry, we add the density structure of the flared disk (equation 3), 
as shown in the top-right panels of Figure 2b and 2c. This disk has a mass of 0.01 M Q and 
an outer radius equal to the envelope centrifugal radius, 100 AU. For our disk parameters, 
the dust destruction radius is 7.5-R* which we set as our inner disk radius. For a disk with 
ddisk = 0.01, and the disk parameters in Table 1, the fraction of luminosity emitted within 
the disk (i.e., beyond 7.5R+) is L acc = 0.0078L*. The third geometry includes the disk and 
a streamline cavity with an opening angle, 9 C , of 25°at 5000 AU. And the fourth geometry 
includes the envelope, disk, and a curved cavity with an opening angle of 20°at 5000 AU. 
Notice Figure 2b shows that the curved cavity carves out more of the inner envelope material, 
even though its opening angle at 5000 AU is smaller than the streamline cavity. Table 1 
summarizes all the model parameters and Table 2 shows the variations for the four models. 

Figure 3 shows a plot of optical extinction vs. polar angle through the envelope to give 
an idea of the optical depths involved. In Model 1, the rotationally-flattened envelope, the 
extinction ranges from Ay = 32 pole-on to Ay ~ 225 at an inclination of 89.999 degrees (the 
solution for density has a singularity in the midplane). In the models with bipolar cavities 
(models 3 and 4), the extinction drops to Ay = 1.5 near pole-on inclinations. In the models 
with disks, the extinctions increase to several million at edge-on inclinations. Thus, the 
star and inner disk radiation will be extincted out to very long wavelengths for the edge-on 
viewing angles in these models. This radiation can escape through scattering at wavelengths 
where the albedo is not negligible. Wood et al. (2002a) shows the varying contributions of 
direct, scattered, and thermal emission for disk models. 



3.1. Temperatures 

Figure 4 shows temperatures calculated from the radiative equilibrium solution in the 
four geometries. Some interesting features to note are: 1) The low-density bipolar cavities 
(models 3 and 4) are hotter than the surrounding envelope, since they are exposed to more 
direct stellar radiation. 4) The midplane is cooler than the envelope or cavity, especially in 
the models with a disk (Models 2-4) (Figure 4b, 4c). Regions in the midplane beyond about 
1 AU see less of the radiation from the star and warm inner regions due to the large optical 
depths through the disk (Figure 3), and thus are heated less than the surrounding envelope. 
4) The interior of the disk in Model 2 is warmer than those in Models 3 and 4 (Figure 4c). 
This model has no bipolar cavity and thus more radiation is trapped within the envelope to 
heat the disk. 

Figure 5 shows the temperature behavior in different regions in more detail. This plots 
temperature as a function of radius for four polar angles, 9 = 0°, 30°, 60°, and 90°. In Model 
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1, which has no disk or cavity, the temperature profiles are similar at all inclinations. The 
equatorial profile is slightly steeper due to the higher optical depths (A v = 225, compared 
to 33 at the poles). The temperature scales roughly as r -05 out to about 3 AU, and then 
transitions to r~ - 33 beyond 200 AU. As shown in Kenyon et al. (1993a), the temperature 
profile in optically thin regions is T oc r~ a , where a = —2/(4 + p), and p is the power law 
exponent for the grain opacity, i.e., k v oc u p . Beyond a radius of 100 AU, the emitting 
temperatures are below 40 K, so the peak of the blackbody function is longward of 75 
fxm. At these wavelengths, the opacity power law exponent is p = 2, giving an optically 
thin temperature profile of r" 1 / 3 in agreement with our results. The light grey line plotted 
with each model is a power law with a = 0.33. In the inner region, we can attempt to 
compare to the diffusion solution, though our opacity exponent varies with wavelength and 
therefore temperature. The diffusion solution (Kenyon et al. 1993a) gives T oc r~ a , where 
a — (1 + n) / (4 +p), and n is the density power law exponent, p oc r~ n . At 1 AU, for the 
rotational infall solution, n = 1/2. The temperature is 400 K, giving peak Planck emission 
at 7.5 fim. At this wavelength there is a kink in the opacity law, where p varies between 
0.5 (on average) to the long wavelength and 1.5 to the short wave. This gives a range of 
a = 0.42 — 0.6 for the calculated temperature profile. Our value of a = 0.5 is within this 
range. Kenyon et al. (1993a) find a value of a ~ 0.7 in similar envelope geometries as our 
Model 1; however, the profiles in Efstathiou & Rowan-Robinson (1991) appear to be closer 
to a ~ 0.5. These differences could be due in part to different grain models and envelope 
parameters. 

The temperature profiles for models with disks and cavities have much more variation 
with polar angle, as shown in Models 2-4 in Figure 5. The midplane temperature profile 
shows a steep drop in the disk, with T oc r~ L15 out to 1 AU, and flattening out to 100 AU. 
Then it follows a similar profile as the rest of the envelope, scaled to a lower temperature. 
The temperature profile in the polar region in the models with cavities (Models 3-4) is similar 
to the models without cavities. It flattens out to the optically thin profile, a = 1/3, at a 
smaller radius, leading to higher overall temperatures in the cavity. The 30°and 60°profiles 
(dotted and dashed lines) diverge from the 0° line where the path intersects the cavity wall, 
and drops to the lower envelope scale. This happens earlier in Model 3 than Model 4 due 
to the shape of the cavity walls. At large radii, the models with cavities are hotter along all 
directions, even in the midplane. 

To better understand the temperature behavior in the disk, in Figure 6 we compare the 
midplane temperature profile of the models with disks (Models 2-4) with those of disks with 
no envelopes. The grey line plotted with each model is a power law, T = 90 (r/lAU)~ a4 K. 
It is clear from these plots that no single power law characterizes the temperature behavior 
in either the disk-only or disk+envelope models. In all of the models except Model 2, the 
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temperatures are similar between about 2-10 AU. In this region, the accretion luminosity is 
not as important, with the disk temperature determined mostly from stellar flux heating the 
flared surface. The midplane temperature is higher in Model 2 than any of the others due to 
extra heating from the surrounding envelope (see also D'Alessio, Calvet, & Hartmann 1997). 
For Models 3 and 4, the cavities are not providing much additional heating to the disk. As 
we will show in the next section, the emissivity in the cavities is actually higher than in 
the same region of the envelope in Model 2. However, this radiation is also able to escape 
through the cavity in Models 3 and 4, whereas the envelope radiation in Model 2 is trapped 
and heats up the disk. Beyond about 10 AU, Models 3 and 4 diverge from the disk-only 
models due to extra heating from the surrounding envelope. In this region, Model 3 is even 
higher than Model 2, probably due to heating of the cavity walls by stellar radiation (due 
to their higher curvature, they intercept more stellar radiation). 

Inside 1-2 AU, the temperature profiles are the same for all the models with disk ac- 
cretion, with a power law exponent a = 1.15. This is higher than expected for passive 
accretion (a = 0.75; Kenyon & Hartmann 1987), but is in agreement with the plotted pro- 
files of D'Alessio et al. (1998) between 0.1 and 2 AU. The temperature is determined by the 
accretion luminosity which has a steep falloff with radius (equation 4). Between 0.073 and 
0.16 AU (7.51 to 16.4 R+), the accretion luminosity raises the disk temperature above the 
dust sublimation temperature of 1600 K. We set the opacity to in this region since the gas 
opacity is much lower than dust. With surface density of 1300-1400 gem" 2 and gas Rosse- 
land mean opacity k r ~ 10~ 4 cm 2 g -1 (Alexander, Johnson, & Rypma 1983), the vertical 
gas optical depth through the disk (top and bottom) is r = ~ 0.1 which is reasonably 
optically thin (Lada & Adams 1992). For the disk without accretion, the temperature profile 
exponent is a = 0.6 from 0.1-1 AU, then flattens out to a = 0.4, and even lower beyond 10 
AU. This is in reasonable agreement with blackbody emission models that give a = 0.5 for 
our flared disk with (3 = 5/4 (Kenyon & Hartman 1987). 



3.2. Emissivity and Infrared Images 

Even though Figures 4 and 5 show that temperatures in the bipolar cavities are higher 
than in the envelope, we might expect the overall emitted energy from the cavities to be 
small due to the low densities. However Figure 7 shows that the high temperatures lead to 
large emissivities in the cavity, despite the lower densities. This is especially apparent at the 
longer wavelengths, 24 and 160 //m, shown here. The emissivity is defined as j u = K v pB u . 
Figure 7 plots azimuthal slices of the emissivity. It shows how much energy per unit volume 
is emitted throughout the envelope (but is not a radiative transfer calculation showing how 
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light is propagated through the envelope). Even though the models with bipolar cavities 
have slightly lower mass than the others, their total emissivities are larger due to the higher 
temperatures in the cavities. Note that the images are plotted to different size scales for 
the three wavelengths. Though the magnitude of the emissivity is much lower at longer 
wavelengths, the brightness profile is much shallower. 

Figure 8 shows images that result from the radiation transfer calculations. These are 
shown at an inclination of % — 80°. From top to bottom, the four model geometries are 
plotted, and from left to right, 3-color composities of near-IR (left), mid-IR (middle), and 
FIR (right) are plotted. The images are displayed to mimimum flux levels described in 
the figure caption, assuming a distance to the source of 140 pc. The near-IR images are 
similar to observations and previous models (Padgett et al. 1999, Stark et al. 2002). The 
mid-IR images should be detectable by the upcoming Space Infarared Telescope Facility 
(SIRTF) in nearby star forming regions, and the FIR images should be detectable by SIRTF 
out to much larger distances. As Figure 7 showed, the mid-IR thermal emission should be 
rather localized near the central source. The large extent of the mid-IR images in Figure 
6 results from scattered light. As discussed in Figure 1, the albedo of the Ry = 4 grains 
used here is higher than ISM dust at these wavelengths. Since the albedo of the dust is 
sensitive to grain size, the mid-IR images should be a useful diagnostic of grain size. The 
FIR images delineate cavity walls and show shadowing in the disk. However, these effects 
will be completely washed out at SIRTF resolution, which, at 70 /iin, has a full width at 
half maximum (FWHM) of 38"or 2380 AU at a distance of 140 pc. We convolved the 
images with Gaussian point spread functions broadened to the HST and SIRTF resolution 
for each wavelength. We don't show them here because the near-IR and mid-IR images are 
essentially unchanged at the resolution shown here, and the FIR images are washed out and 
indistinguishable from each other, though they are elongated with about a 2/1 ratio. 



3.3. Spectral Energy Distributions 

The model SEDs are shown in Figure 9. Ten inclinations are plotted, the top curve 
corresponding to pole-on and the bottom to edge-on. The absorption feature at 3.1 fim is 
due to ice opacity, and that at 10 /zm is due to silicate. A distinctive feature of the SED plots 
is the broad dip at 10 /im at low inclinations in models that include the flared disk. This is 
due to the disk geometry and the grain properties: At shorter wavelengths, even though the 
opacity through the disk and envelope is larger (Figure 3), the albedo increases allowing light 
to scatter out at all inclinations. Note that the near-IR colors are determined completely 
by scattering for inclinations i > 60 (cosi < 0.5), so the concept of "de-reddening" based 
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on the colors is meaningless since that only applies to direct extincted light, not scattered 
(see Stark et al. 2002). At wavelengths longward of 10 /im, the opacity decreases enough to 
allow thermal emission to propagate through the disk and envelope. Near 10 /im, the low 
albedo and high opacity combine to prevent radiation from escaping at high inclinations. 

Another striking feature of the SEDs in Figure 9 is that the models with bipolar cavities 
(Models 3 and 4) are much brighter at optical/near-IR wavelengths in better agreement with 
observations than the models without cavities (Kenyon et al. 1993a,b). Notice in Models 3 
and 4 that the pole-on inclinations still suffer extinction at short wavelengths due to the dust 
in the cavity. The shortwave clearly is very sensitive to the amount of dust in the cavity. 
Models 3 and 4 also show the 10/xm silicate feature in emission at pole-on inclinations, due 
to low optical depth to the source of emission (the inner disk surface). 

It is clear from the plots in Figure 9 that the more edge-on spectra have lower integrated 
fluxes than other inclinations. In a spherical source, the integral of the measured SED 
multiplied by And 2 , where d is the distance to the source, is equal to the luminosity of the 
source. However, this is not true for sources with non-spherical envelopes, where the SED 
must be integrated over direction as well as frequency, i.e., 



where \i = cost and i is the inclination, also the polar angle. We note that in our models, this 
integral is equal to our input luminosity to within 0.02% (the difference is due to performing 
the integral after coarsely binning the results in angle, since flux is fully conserved in our 
models). Figure 10 shows the correction factor that would have to be applied to the inte- 
grated SED at a single inclination to get the true luminosity of the source. The correction 
factor is simply the luminosity divided by the integrated SED at each inclination i, 



As Figure 10 shows, the observed SEDs can estimate the luminosity of a source only to 
within a factor of 2. While it is not surprising that the edge-on integrated flux is lower than 
the true luminosity, it is interesting to note that the pole-on integrated flux is a factor of 2 
larger than the true luminosity for the sources with bipolar cavities. This inclination receives 
a maximum of direct stellar flux, scattered flux and reprocessed emission. 

Nakazato et al. (2003) defined a parameter, f L) as the ratio of the integrated SED, 
J F v dv to the peak flux uF u , and argued that this can be used as an indicator of inclination. 
We had to modify their definition slightly because for our models with cavities, the wave- 
length of the peak shifts from about 70 /im to 1.6 /mi at the pole-on inclinations, causing 




(7) 



correction = 



lF v {i)dv 



(8) 
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fi to decrease for pole-on inclinations. We can resolve this by defining f L (FIR) to be the 
ratio of the SED to the peak flux in the far infrared (FIR). Figure 11 plots f^FIR) vs. 
inclination for our set of models. Our models also show a correlation with inclination, but 
the behavior is different from Nakazato et al.'s values. This is due to the differences in our 
model geometries. Nakazato et al. include an envelope, disk and cavity in their models but 
their envelope is spherically symmetric up to the cavity opening angle. Thus, their fi is 
relatively flat up to the cavity opening angle and then flat within the cavity at a high value. 
In contrast, our models using a rotationally flattened envelope show a more gradual increase 
in f^FIR) as the inclination becomes more pole-on (increasing cosi), due to varying en- 
velope extinction with polar angle, a different cavity shape, and a different assumption for 
the dust content in the cavity. Thus, the behavior of Jl with inclination depends strongly 
on model assumptions. Nakazato et al. claim this parameter fi can determine inclination 
to ±5°, and use as an example IRAS 04365+2535 for which they find an inclination of 20°. 
However, the high near-IR polarization (Pk ~ 10 — 15%; Whitney et al. 1997, Lucas & 
Roche 1998), and molecular outflow observations (Chandler et al. 1996) suggest inclinations 
of i ~ 40° — 70°. We suggest that more detailed modeling, including as much observational 
information as possible, is required to "calibrate" the models before applying this factor 
based on SED information alone. 



3.4. Color-Color Plots 

Figure 12 shows color-color plots for our models. The black crosses are main se- 
quence/giant branch stars, and the dotted lines extending from them show the reddening 
lines. The region enclosed by these lines shows where the MS/GB falls on these plots. The 
model points are plotted with different shades of grey representing inclination, with light 
grey representing pole-on, and black edge-on. The models without cavities are shown by 
the filled circles (Model 1) and asterisks (Model 2). The models with cavities are shown 
by open triangles (Model 3) and open diamonds (Model 4). The models with cavities are 
substantially bluer in all the plots. In a source observed to have a bipolar cavity (as indi- 
cated by images resembling Figure 8, outflows, and high polarization), it is clear from Figure 
12 that 2-D cavity models should be used to interpret their colors. Using 1-D models, or 
models without bipolar cavities, would lead to assignment of a later evolutionary state to 
the source than is appropriate. Another interesting feature in Figure 12 is that some of the 
edge-on sources are the bluest. This is because the disk blocks all the direct light, leaving 
only scattered light which is blue, especially in the MIR where the albedo increases steeply 
towards shorter wavelengths. 
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The bottom left panel of Figure 12 shows that the Class I source colors are well-separated 
from reddened main sequence and red giant stars for this particular color choice in the mid- 
IR and our grain model. The bottom right panel, plotting [3.6]-[5.8] vs [8.0]-[24.0] shows 
that Class I sources are very red for this color choice. The main sequence stars have to 
be reddened by at least Ay > 100 to overlap with these sources. Thus, in nearby star 
forming clouds, where foreground main sequence stars have almost no reddening, the Class 
I sources will likely stand out. Our next paper will explore in more detail how well different 
evolutionary states are separated in color-color space. 



3.5. Polarization 

An additional tool to distinguish the various geometries as well as grain properties is 
the polarization spectrum, shown in Figure 13. Plotted is Q = Pcos(x), where x is the 
orientation of the polarization vector. Because these models are axisymmetric, and the 
grains are spherical (no dichroism), the integrated or net polarization is either parallel or 
perpendicular to the rotation axis. If Q/I is less than 0, the net polarization is perpendicular 
to the rotation axis. This is the case over much of the spectrum because most of the scattered 
light emerges from the bipolar cavity regions. In all of the models, the polarization is highest 
for the edge-on inclinations and decreases to pole-on. The polarization is sensitive to the 
envelope optical depth: To get high polarization from scattered light requires optical depths 
to the unpolarized stellar/disk radiation of r > 5 (Whitney et al. 1997). We see this 
transition in the 1-10 fim region, depending on inclination. At longer wavelengths, more 
scattered radiation emerges from the equatorial region because the optical depths in the 
polar regions are too low. We should expect to see the polarization reverse sign at these 
wavelengths. However, the grain albedo is decreasing rapidly in this region so the polarization 
vanishes before reversing sign (though there is a hint in Model 4 of a polarization reversal). 
In a less embedded source we would expect the polarization to reverse sign at a shorter 
wavelength so the polarization reversal would occur before the albedo drops to 0. We will 
explore this in more detail in future papers. In disks, where the grains are larger, the albedo 
should remain higher towards longer wavelengths (Wood et al. 2002a). Thus, polarization 
spectra should be a useful diagnostic of grain growth as well as evolutionary state. 

Figure 13 also shows the behavior of polarization across the ice band feature at 3 /im. 
The polarization increases across these features because the opacity increases. With higher 
opacity, less direct (unpolarized) light is present to dilute the polarized light. This behavior 
agrees with observations across the 3/xm feature (Hough et al. 1996, Kobayashi et al. 1999, 
Holloway et al. 2002). 
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4. Conclusions 

We have shown that 2-D models of realistic Class I source geometries produce substan- 
tially different SEDs from 1-D and simplified 2-D models. The inclusion of a flared disk 
and bipolar cavity changes the SEDs and colors substantially. The effect of the disk is to 
cause a broad dip near 10 jum in nearly edge-on sources (i > 80°), and the cavity allows 
more shortwave radiation to escape at all inclinations. Both the disk and cavity introduce 
inclination effects much greater than the simple rotational infall model gives (Model 1). In 
addition to the spectral shape, the integral of the SED also varies with inclination. In 1-D 
models, this integral is equal to the luminosity of the source. In 2-D or 3-D models, the 
flux must be integrated over direction as well, as shown in equation (7). This is impossible 
to do for an observed source which is only viewed from one direction. Figure 10 shows the 
correction that would have to be applied to a source with the same geometry as our models 
to get the true luminosity of the source for a given inclination. 

Nakazato et al. (2003) claim that the inclination can be determined from a factor fi 
which is the ratio of the integral of the SED to the peak of F v dv. Our models give a different 
behavior for fi due to the different envelope geometry. We conclude that fi on its own 
is not a robust indicator of inclination. Perhaps incorporating many more observational 
constraints to the models could provide some calibration to this factor. 

We placed our model results on color-color plots (Figure 12) for easier comparison to 
observations, and find that the models with bipolar cavities are much bluer than those 
without. In addition, the models with cavities show a large spread with inclination. It is 
clear from these plots that realistic 2-D models are required to properly interpret the colors 
of Class I sources with bipolar cavities. A 1-D model would require a lower envelope mass 
to fit the same colors, which would imply a later evolutionary state for the source than is 
appropriate. 

The images in Figure 8 and polarization spectra in Figure 13 can be used to understand 
dust properties in star forming regions. Our dust model, fit to the extinction curve typical in 
the Taurus molecular cloud, has enough scattering albedo in the mid-IR to produce extended 
scattered-light images in the mid-IR. As shown in Figure 1, this albedo is higher than the 
ISM value. These images are most easily detected at edge-on inclinations where the flared 
disk blocks the bright direct stellar and inner-disk radiation. The sensitivity of SIRTF may 
be enough to image these sources and confirm the higher albedo due to larger grains. The 
polarization spectra also provide constraints on both envelope geometry and grain properties. 
Our choice of an ice mantle occupying 5% of the outer radius of the grains gives a 3.1 //m 
feature in good agreement with observations. Further models at different evolutionary states 
will be needed to better characterize the polarization behavior in protostars. 
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Given that bipolar cavities and inclination effects introduce a large spread in observa- 
tional parameters for Class I sources, we suggest that images, outflows, and polarization 
measurements can aid in determining cavity geometry and inclination. Then 2-D models 
can more accurately assign envelope mass and evolutionary states to Class I sources. In our 
next paper we will explore in more detail the effects of evolution on the SEDs, colors, images, 
and polarization of protostars. 

We thank Martin Cohen for providing the colors of main sequence/giant branch stars 
in Fig. 12, as well as providing zero-point magnitude scales for several SIRTF wavebands. 
This work was supported by the National Science Foundation (AST-9909966, AST-9819928) 
and the NASA Astrophysics Theory Program (NAG5-8587). K. Wood acknowledges support 
from the UK PPARC Advanced Fellowship. 
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Table 1. Model Parameters 



Parameter 


Description 


Value 




ltar ramiic 
O Lclldl 1 dUlUo 


9 nq R 


T 




zinnn v 


T 


Stellar luminosity 


i r 

1 Lq 


M 


oieiiai mass 




1\/T 

J- Vl disk 


Disk mass 


U.Ul lviQ 


a 


Disk radial density exponent 


2.25 


(3 


Disk scale height exponent 


1.25 


h 


Disk scale height at R* 


0.01 


Rd,m'm 


Disk inner radius 


7.5 R* 


Rd,max 


Disk outer radius 


100 AU 


ttdisk 


Disk Viscosity parameter 


0.01 


Lace (Rd,min) 


Disk accretion luminosity to Rd,min 


0.0078 U 


Mdisk 


Disk accretion rate 


1.4 x 10- 8 M Q /yr 


Re,min 


Envelope inner radius 


7.5 -R* 


-Re, max 


Envelope outer radius 


5000 AU 


M c „v 


Envelope infall rate 


5 x 10~ 6 M /yr 



Table 2. Model Variations 



Name 


Disk 


Cavity 


Cav. shape 


0c 


Env. Mass 


Model 1 


n 


n 






O.189M 


Model 2 


y 


n 






O.189M 


Model 3 


y 


y 


Streamline 


25° 


O.177M 


Model 4 


y 


y 


Curved 


20° 


O.176M 
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Fig. 1. — Behavior of the Ry = 4 grains with wavelength. For comparison, the dotted lines 
show behavior of ISM grains (Ry = 3.1). Top left: opacity, k; top right: albedo, u; bottom 
left: average cosine of scattering angle, g; bottom right: maximum polarization, p ; . 
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Fig. 2. — Density plots showing the four geometries. Model 1: rotationally flattened infalling 
envelope; Model 2: same as Model 1 + flared disk; Model 3: same as Model 2 + streamline 
cavity; Model 4: same as Model 2 + curved cavity. The scaling is logarithmic, plotted to 
six orders of magnitude from the peak shown. The contours each decrease by 1 order of 
magnitude, and match the labels of the color bar. The plotted regions and peak densities 
are (a) 5000 AU, 3 x 10" 14 gm cm" 3 , (b) 300 AU, 1.5 x 10" 12 gm cm" 3 , and (c) 20 AU, 
8 x 10~ 12 gm cm~ 3 . 
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Fig. 3. — Optical extinction, Ay, from the central star through the envelope vs. inclination 
with respect to the polar axis, i. The dotted line shows extinction due to the envelope, the 
dashed line is that of the disk, and the solid is the sum. Note the large extinctions at edge-on 
inclinations for models with disks (Models 2-4). 
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Fig. 4a. — Temperature plots for the 4 geometries. Images are scaled logarithmically with 
the color scale shown at the right. The contours correspond to the labelled values on the 
color bar (in units of log(T)). The plotted regions and temperature scales are (a) 5000 AU 
and T = 8 - 300°K, (b) 300 AU, T = 22- 1000°K, and (c) 20 AU, T = 42 - 1600°K. 
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Model 1 Model 2 




0.1 1.0 10.0 100.0 1000.0 0.1 1.0 10.0 100.0 1000.0 

R(AU) R(AU) 

Fig. 5. — Temperature as a function of radius along four co-latitudes, 9, with respect to the 
polar axis: 9 = (dot-dashed line), 9 = 30 (dashed), 9 = 60 (dotted), and 9 = 90 (solid). 
The light grey line is a power law, T = 180 (r/lAU)~ a33 K 
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Fig. 6. — Temperature as a function of radius in the midplane of selected models with disks. 
Each model below the top is scaled down by a factor of four from the one above. The grey 
line is a power law, T = 90 (r/lAU) -0 ' 4 K, also scaled down by a factor of four for each model. 
Notice that the disk midplane is warmer throughout the Model 2 disk due to radiation from 
the surrounding envelope. 
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Fig. 7. — Emissivity plots at three wavelengths: 8 /im (left), 24 /im (middle) and 70 jum 
(right). The four model geometries are plotted from top to bottom. The scaling is logarithmic 
with the color bar shown on the right, scaled to the peak. Peak emissivities plotted are 
1.5 x 10 _1 at 8 /im, 1.5 x 10~ 4 at 24 //m, and 3 x 10 -9 at 160 /xm, in units of erg/cm 3 /s/Hz/sr. 
From left to right the region plotted spans a radius of 20, 300, and 5000 AU. The emission is 
more centrally concentrated at shorter wavelengths and in models without cavities (Models 
1, 2). The contours each decrease by 1 order of magnitude, starting from the peak. 
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Fig. 8. — Three-color composite images of the four models (going from top to bottom) in 
three different wavelength regions (left to right), viewed at an inclination i = 80°. The near- 
IR images (left) are composites of images calculated through three NICMOS filters centered 
at 1.1, 1.6 and 2.05 fim. The mid-IR images (middle) are composites at three SIRTF IRAC 
filters centered at 3.6, 4.5, and 5.8 /im. And the far-IR images (right) are for the SIRTF 
MIPS filters centered at 24, 70 and 160 /im. The images are scaled logarithmically with 
minimum intensities of 0.01, 0.0074, and 0.5 MJy/sr for the near-IR, mid-IR, and far-IR 
images, respectively, assuming a distance of 140 pc. 
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Fig. 9. — Spectral energy distributions for the four geometries. Fluxes (erg cm -2 s _1 ) are 
scaled to a 1 L Q source at a distance of 140 pc. Ten inclinations are plotted ranging from 
cosi = 0.05 to 0.95 at intervals of 0.1. The light grey line is the input stellar spectrum. 



-34- 



u 
o 



c 

o 

o 

CD 

o 

CJ 



2.5 



2.0 




0.5 



0.0 



0.0 



0.2 



0.4 0.6 
cos(i) 



0.8 



1.0 



Fig. 10. — Luminosity correction factors as a function of inclination for the four geometries: 
Model 1: dot-dashed; Model 2: dashed; Model 3: dotted; Model 4: solid. These curves 
show the correction factor that would have to be applied to the integrated SED at a given 
inclination to correctly determine the source luminosity. 



-35 - 




Fig. 11. — The ratio of the integrated SED to the maximum far-IR flux is plotted against 
inclination. Model 1: dot-dashed; Model 2: dashed; Model 3: dotted; Model 4: solid. 
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Fig. 12. — Color-color plots. The symbols correspond to the four models as follows, Model 
1: filled circles; Model 2: asterisks; Model 3: open triangles; Model 4: open diamonds. The 
models with cavities (3 & 4) are much bluer than the models without (1 & 2), despite their 
similar envelope masses (Table 2). The symbols are plotted at different levels of greyscale 
to indicate inclination, as shown by the greyscale bar in the top left panel. Edge-on sources 
are black, and pole-on are light grey. Notice that some of the bluest sources are edge-on. 
This is due to the dominant contribution of scattered light to the source flux. These edge-on 
sources are very faint so would stand out in a color- magnitude diagram as faint and blue. 
The dotted lines are reddening vectors off the main sequence and red giant branch (small 
crosses). The black dots along one of the reddening lines denote increments of Ay = 30. The 
models in the bottom left panel are well-separated from reddened main sequence and giant 
stars. The models at bottom right are well-separated from stars with Ay less than 100. 
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Fig. 13. — Polarization spectra for the four geometries. Five inclinations are plotted from 
cos 2 = 0.05 to 0.85 in intervals of 0.2. The lines are lighter at lower inclinations. 



